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Abstract 

The methods commonly used for numerical differentiation, such as the "center-difference formula" 
and "four-points formula" are unusable in simulations or real-time data analysis because they require 
knowledge of the future. In BardTl, an algorithm was shown that generates formulas that require 
knowledge only of the past and present values of f(t) to estimate fit)- Furthermore, the algorithm can 
handle irregularly spaced data and higher-order derivatives. That work did not include a rigorous proof 
of correctness nor the error bounds. In this paper, the correctness and error bounds of that algorithm 
are proven, explicit forms are given for the coefficients, and several interesting corollaries are proven. 
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1 Introduction 

In numerical analysis, one is often taught that for positive h, rather than use 

f(t) - f(t - h) 



it is better to use 



/'(*) 



/'(*) 



f(t + h)-f(t-h) 



2h 



(the "backward difference" formula) 



(the "center difference" formula) 



or better still 



f( t +h)-f(t-h) 

2h 



f(t + 2h)-f(t-2h) 
Ah 



(the "four points" formula) 



because of the error of the approximations. The error in the first case is 0(h), compared to 0(h 2 ) in the 
second case and 0(h 4 ), where 0(h m ) indicates some function that is bounded when divided by h m for all 
sufficiently small values of the step h. 

In simulations or in real-time data analysis, these more advanced formulas are of no use, because they 
require knowledge of the future — i.e. f(t + h) and f(t + 2h). However, the future is either unknown (in the 
case of real-time data analyses) or not yet calculated (in the case of simulations). 

In a previous paper [1], Bard showed how to compute similar formulas such as 



)f (t - 4h) -tf(t- 3h) + 3/ (t - 2h) -4f(t-h) + (t) = f (t) ^ 4 / (5) (*) + Wf^ (t) 



4 , \ / 3, v / / , v / 12 , v , , v , 5 , w 3 

that work with only the knowledge of the past and present, but that do not require knowledge of the future. 
Furthermore, the algorithm presented there also allows for irregularly spaced sampling of /(t), and arbitrary 
numbers of data points, as well as second and higher derivatives. The error bounds of the formulas from that 
paper were handled heuristically, and only a sketch of the proof of correctness was given. Here we give a 
rigorous proof, explicit forms of the coefficients of those formulas, and several rigorous statements about the 
error bounds. The obtained formulas could be used to reduce the running time while keeping the accuracy 
of derivative estimates. 

For a survey of previous work in this topic, see Section IX of [1]; that paper also contains several extended 
examples, and SAGE code for carrying out the algorithm. This document can be thought of as a sequel to 
that paper. 
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2 Initial Calculations 

Assuming that / is at least (n — 1) times differentiable, we use the following form of Taylor's Theorem 



/(*+*)= £/ (i) (4 



R„ 



where represents the ith derivative of /(i), and i?„ is the remainder term in one of the well-known forms. 
Two of those forms will be explored in detail during Section 4. We define the Oth derivative of a function to 
be the function itself. 

Let us suppose that 81,62, ■■■ ,8 n are distinct real numbers. Using the formula above we arrive at the 
following representations for f(t + 8jh) 



8% l 



(1) 



For any integer k such that < k < n, any real non-zero h we are looking for a linear combination of 
terms f(t + Sjh) that approximates f^ k \t) term up to some reasonable remainder term. More precisely, we 
start with the linear combination Ci/(t + Sih) + ■ ■ ■ + c n f(t + S n h), rewrite it using (JTJ, and look for the 
values of Cj that eliminate as many terms as possible except for the terms involving f^ k \t). 



5>/(t + *i/0 = E 



3=1 



3 = 1 
n-1 



v i=0 ' / 



i=0 



J=l 



c 3 Rn 

3 = 1 



(2) 



(3) 



It is useful to note that in the above we interchanged two summation symbols, which was legal because 
both sums are finite; we also pulled out of the summation over j any factors that do not depend on j. So 
far, everything we have said is true for all real CjS. 

If the numbers ci, . . . ,c n are chosen so that the following conditions are satisfied 



3=1 

n 

3=1 

then ([3]) dramatically simplifies to become 



for < i < n — 1 , i ^ k 



3=1 3=1 



which yields the desired result 



-. 71 -. 71 

y; E c ^ + = / (fe) (*) + T* E c ^ 

3=1 3=1 

Accordingly, to simplify matters, we set 

^ n ^ n 

R = —f^Y. c 3 R ^3 therefore ^ ^ Cj f(t + 8 s h) = /<*>(*) + R 



(4) 
(5) 

(6) 

(7) 

(8) 



3=1 3=1 
Thus R represents the error of the formula, and we will seek to place bounds on R in Section |U 



2.1 Existence and Uniqueness 



We still have to prove the existence of numbers ci, . . . , c n that satisfy ^ and ©. In fact, we will prove that 
such numbers exist, and that they are unique. The equations represented by Q and ([5]) could be rewritten 
as a system of linear equations Ac — b with = (5* _1 and = fc! with the rest of coordinates of b 
are zero. Matrix A is an example of a Vandermonde matrix; A is invertible precisely if and only if all of 
the 8i,...,S n are distinct (see Equation (|lip ). Because matrix A is invertible, there exists unique solution 
c = A~ 1 b, which is exactly what is required. In fact, vector c is a column of matrix A^ 1 multiplied by fc!. 
The above calculations could be stated as the following algorithm. 

2.2 The Algorithm 

Input: Any n distinct real numbers Si, 82, ■ ■ • , 8 n , and an integer fc such that < fc < n. 

Output: A formula for the fcth derivative of any n-times differentiable f(t), in terms of f(t + 8jh), for any 
positive real h, with error proportional to h n ~ k and the nth derivative of f(i). 

1. Define the n x n matrix A such that Aij = • 

2. Define the n-dimensional vector b to be all zeros, except the bk+i = fc!. 

3. Let the n-dimensional vector c be the solution to Ac = b. 

4. Return the formula 

3=1 

2.3 The Condition Number of A 

We would like to note that because A is a Vandermonde matrix, some readers maybe concerned about the 
numerical stability of solving Ac = b. This is because the condition number of a Vandermonde matrix can 
become extremely tiny as any two data points move close together. However, this is not an issue, because 
in any practical situation numbers Sj would be integers that do not depend on the value of the step h. 
Furthermore, the value of n would satisfy n < 20. Moreover, using computer algebra packages, one can solve 
Ac = b using exact rational arithmetic in less than a second, totally avoiding any floating-point computations 
of any kind. Therefore, the condition number of A is immaterial, and c will be known exactly. 



3 Explicit Formulas for the qs 

While the algorithm will very quickly produce the correct values of c, it might be useful to have some explicit 
formulas in order to deduce properties of the CjS. We will accomplish this by first finding some formulas 
for the determinants of the minors of a Vandermonde, and then use those with Cramer's rule to obtain the 
explicit formulas for the qs. 

3.1 Some Useful Notation 

We define a particular n x n-matrix V(yi, y%, . . . , y n ) via Vy = Vj~ ■ 

Recall that A = V(Si, S2, ■ • • , 5 n ) in our algorithm. We will abbreviate this A = V(S). Furthermore, we 
will use a hat to indicate the removal of one of the deltas. Explicitly, 

Sj = (81, . . . , Sj-i,Sj + i, ...,S n ) 

Then it shall be useful to denote by A the largest absolute value of any 5. Let ej be the distance from 
Sj to the nearest other S. Finally, let e signify the smallest of all the CjS. Explicitly, 

A = max{|di|, . . . , \S n \} and Cj = min \5j — 8i\ with e = min{ei, . . . , e n } (9) 



Since we assumed that h > 0, and 5%, . . . , S n are real numbers, we define the interval of consideration I 
to be the smallest closed interval containing all of (t + Sih) values: 

I = {x : min (t + Sih) < x < max (t + Sih)} (10) 



3.2 The Determinants of the Minors of a Vandermonde Matrix 

In any case, the determinant of V, which we denote v, has a well-known formula 

v{yi,V2, ■ ■ ■ ,Vn) = n (Vj-Vi) (H) 
l<i<j<n 

which is simply the product of the differences of all possible distinct pairs of yi and yj , taking care to always 
subtract the y of higher index from the y of lower index. 

The following lemma establishes a formula for the determinant of a minor of a Vandermonde matrix. We 
presume that this must have been known for quite some time, but we could not find a proof of it anywhere, 
and we find the following proof both short and simple. 

Lemma 1 Consider V{y\,y2, ■ ■ ■ ,yn), with all ys distinct. Let M{i,j) denote the minor formed by deleting 
the ith row and jth column from V . The determinant of M(i,j) is given by 

m(i,j) = det M(i,j) = v(yj)<T n -i,n-i(Vj) (12) 
where o n -\. n -i is the (n — i)-th-degree symmetric polynomial in n — 1 variables 

The proofs of the lemma can be found in Appendix [XJ of the full version of the paper on arXiv . org. 

3.3 Cramer's Rule 

Now suppose we want an explicit formula for Cj in c, as produced by our algorithm. As in our algorithm 
and Lemma 1, we define A = V(5) . Using Cramer's Rule, we can take A but replace column j with 6; the 
determinant of that modified matrix, divided by the determinant of the original A, equals the value of Cj. 
We should expand the determinant of the modified matrix on column j, which is equal to b, because b is 
zero in all but one entry. We would then obtain: 



{-if+iAtj detMy + (-l) 2 +^A 2j detM 2j + ■■■ + (-l) n+j A nj detM n 



det A 

v(S) 



= ( lf+i+l h A 5 3) (J n-l-,n-k-l{0~ j ) 



However, we should observe that 

v(6j) _ 



v & (rw^,-^)) (ii r , <\. 



because all terms in the denominator not involving Sj will also be found in the numerator. Therefore, we 
can conclude ^ 

c . = r_i\k+j+i k\a n - ltn - k -i(Sj) 

(II: , ,^ <>< : )0L< ,A»> 

or simply 

{-l) n - k kla n -i, n - k -i{5j) 



11^ is. Si) 



When we specialize the above formula to the case Si = —i, for k = 1 we obtain 
and for k = 2 we obtain 

The proof of the special cases can be found in Appendix [Ej of the full version of the paper on arXiv . org. 

4 Error Bound Theorems 

In this section we return to the formula used in the main algorithm, and make the error terms more explicit, 
and provide some useful error term estimates. 

4.1 A Simpler Form of a Particular Sum 

We are about to use the sum S = X^=i c i^j m the following discussion, so we start by working on a simpler 
closed form of the sum. 

Lemma 2 The sum S has the following closed form: S — Yl^—x Cj5 r f = (— l) fc+n+1 fc!(T„. n _fc((5) 
We will need the following estimate on S as well 

Lemma 3 The sum S — V™ , c,-<5? satisfies \S\ < > |c,<5™ < =- — ^7 

j=l 

The proofs of the lemmas can be found in Appendices |B] and [C] of the full version of the paper on 
arXiv.org. 

4.2 Error Term Estimates 

The easiest to remember form of the error terms in (JTJ is the Lagrange form of the error term: 

Rn, = ™> (14) 
for some £j in the interval between t and (t + Sjh). Using equation we immediately have 

1 " f( n Mf )5 n h n h n ~ k n 

R =-^J2 * „, 3 = -— E (15) 

4=1 ' ' 3=1 

While this form of the error term is correct and exact, it is hard to use in practice, because the £s are 
unknown. 

Theorem 1 If 

1. f is n-times differentiable on an open set containing interval of consideration I , 
2- |/ (n) (r)l < M f or all r G I, with I defined in (TO]) 
3. A, Cj, e as defined in ([9]) 
then the algorithm's error term (see Equation ([SJj satisfies the following estimate 

l/A2n-/:-l 

\R\ < -h n (16) 

' ' ~ £"-!(«- k- 1)! V ' 



Proof. This is a direct calculation starting with (fT5|) . and using Lemma [31 

i n - k 

n\ e n - l (n-k-l) 



ijn—k n Kn — h 71 un—k A 2n — k — l rj | 

1^1 < ^Y\cj$U in) (M < M^—Y\c^\ < M-— % 7 (17) 

3=1 j=i 



Rearranging terms we arrive at the estimate we're looking for. □ 
We note that another easy to remember form of error term in Taylor theorem is the little-oh form. In 
our context, we rewrite (JJJ) in the form 

f (n Ut)S r 'h n 

Rn,i = i + oa fl (18) 

We must confess that in the above we wrote an extra term in the Taylor expansion for future use plus the 
actual error term. Since only h is typically considered as a parameter approaching zero, and Sj is a constant 
we can replace o(h n 5™) with o(h), and using that, we write (JSJ as 

R = ~hk2^ c i y n\ + ( } J = n\ j j ) 

^-^""^""^'r^—^+o^-^) (19) 

where we used Lemma [5] and the o(h n ~ k ) term is the sum of n different o(h n ~ k ) terms coming from (|18l) . 
We state this result as the following theorem. 

Theorem 2 If f is n-times differentiable on an open set containing interval of consideration I then the 
algorithm's error term in (JSJ is 

i ,^(-i)"- fc ^- fc ^/ ( " ) w^.- fc (fl +0( ,„- fc) (20) 

nl 

Sometimes it is useful to know if the algorithm's formula will always underestimate or always overestimate 
the value of f^ k) (t). Using Theorem [2] we are able to prove the following additional result. 

Theorem 3 Assume that all Sj < 

1. If f (n) (t) > 0, then the error term R is positive for all sufficiently small values of h, and the estimate 
of f^ k \t) provided by ([7]) is an underestimate. 

2. If f {n) {t) < 0, then the error term R is negative for all sufficiently small values of h and the estimate 
of f k '(t) provided by ([7]) is an overestimate. 

The proof of the theorem could be found in Appendix [Dl of the full version of the paper on arXiv . org. 
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A Proof of Lemma 1 



Repeat of Lemma 1 Consider V(y\, y%, y n ), with allys distinct. Let M(i,j) denote the minor formed 
by deleting the ith row and jth column from V. The determinant of M(i,j) is given by 



m(i,j) = det M(i,j) = v(yj)a n -i !n -i(yj) 
where o n -\, n ~i is the [n — i)-th-degree symmetric polynomial in n — 1 variables 



(21) 



Proof. Instead of focusing on the single formula for the value of m(i,j) we will find all values of 
m(l, j), . . . ,m(n, j) in one calculation. We start by setting up the matrix V(y±, . . . , yj-i, x, j/j+i, . . . , y n ) 



V(yi,---,Vj-i,x,yj 



,Vn) 



1 

Vl 

III 



i 



i 



Vj-i ■'■ 



1 

II 2 

Vj+i 



1 



x Vj+i 

2 „,2 



T n-l ,.n-l 



y 



n-l 



(22) 



and observing that the determinant of the matrix could be obtained via expansion along the jth column to 
yield an (n — l)th degree polynomial in x: 

Pn-i(x) = v(yi, . . .,yj-i,x,yj+i,. ..,y n ) = a + a%x H h a„_ix ,l_1 = 

(-l) 1+ 'm(l, j) + (-l) 2+ %i(2, j)x + ■■■ + (-l) fc+ %i(fc, j)x h ~ l + ■■■ + (-l) n+j m(n,j)x n - 1 (23) 



From the above, we have an expression for the coefficients of the polynomial: 

a fc _! = (-l) k+j m(k,j) 



(24) 



Note that if x = yu for k ^ j then columns k and j are identical, the determinant is zero, and therefore 
p n -i(x) has a root at x = yk- Because p n _i(x) has degree (n — 1), and all the ys are distinct, this is a 
complete list of the roots: 



x = y n 



(25) 



x = yi, x = y 2 , X = yj-x, x = y j+1 , . 

we immediately have the following alternative expression for p n _i(x) 

Pn-i(x) = a n -i{x - yi) ■ ■ ■ (x - yj-i)(x - y j+1 ) ■■■(x-y n ) (26) 
removing all the parentheses we obtain yet another form of the same polynomial, made from Vieta's formulas: 

Pn-l(x) = a„_iX ,l_1 + a„_i(-l) V n _ M (y.,)ir"~ 2 H 

+ a n _i(-l)"- fe CT„_ 1 , n _ fc (y j )^" 1 + • • • + On-i(-ir~V n _i, n _i(tfj) (27) 
So, we arrive at still another expression for the coefficients of the polynomial: 

<Zfc_i = a n -i(—l) n ~ o- n -l,n-k{yj) 
Equating the two expressions for the coefficients, we immediately obtain a formula for m(k,j) 

m(k,j) = (-1)*+Ja fc _ x - (-l) fe+J a n _ 1 (-l)"- fc (T„_ 1 , n _ fe (%) - (-l)"^a n _ 1< 7 n _ 1> „_ fc (%) 
since by we know 

a„-i = (-1)" +J m(n,i) = (-l)"+^(%) 

we have therefore 

m(k,j) = (-l)" +;7 a n _i(T n _i,„_fc(yj) = (-l)" +J (-l)™ + ^(yj)CT n _i !n _fc(?/j) = u(yj) CT n-i,n-fc(j/j) 



(28) 
(29) 



as we claimed. 



(30) 
□ 



B Proof of Lemma 2 

Repeat of Lemma 2 The sum S has the following closed form: 

n 

S = Y. °i 5 l = {-i) k+n+1 ky n ,n- k (s) (31) 

Proof. Using the formula for Cj in its not yet fully simplified form from (|13[) . we see that 

s = i>a? =f:^(- i )^ +ifc! ^ K y" i(gj) - ^E^(-i) fc+j+i -(^)-»-i,n- fe -i^) (32) 

at this point we observe that 

SjO'n-l,n-k-l(Sj) = 0- n , n -k(5) - <7 n _i,„_fc((5j) (33) 

as a general property of symmetric polynomials so that we can continue with simplification of S: 
, l n 

S = -^E^ 1 ^ 1 )^ 1 ^^ r".-^ 1 ) - ^-1,-^(1)) ( 34 ) 
v \ d ) j= i 

«(*) ^ "(*) ^ 

Now, we observe that the first sum in (|35[) is almost an expansion of a Vandermonde determinant along the 
last row, since 

n 
i=i 

n 

= Y^i-^srMs^n-iAl) (37) 

n 

= ^c-ir+^rv?,) (38) 

i=i 

(39) 

which gives us 

^^'(-il^K^W = ^-i, n -fc(*)^*i(-i)* H - 1 - B E(- 1 ) n+it '^)^- 1 (40) 

= ^.^^-L^-l)^-"^) (41) 
v{6) 

= (-l) fe+1 -" ( r n _ lin _ fc (5)fc! (42) 

Having dispatched the first sum, we can now consider the second sum from (|35|) : it involves the coefficients 
for Cj for the (k — l)-th derivative. We use the original formula for Cj for the k-th derivative from (fT3"|) which 
is 



c 



(fc) 



1 



(-ljfc+i+i*!,,^.)^! fc.^tf.) (43) 



to obtain 



cf- 1] = (-l) k+ i{k yiv&a^n-tfa)^ (44) 



and observe that the second sum in (|35[) can become 



7 l n n 
— ^-^-lf+i+^Vn-lin-kfa) = -kJ^(k-l)\S«-\-l) k+ iv(%)* n ^ n - k (%)— (45) 



-ijfc-i) 



-fc^ An, 
3=1 

-kbl k - r) = 



(fc-i) 



(46) 

(47) 
(48) 



for all values of k that satisfy < k < n. We used b n to denote the last coordinate of the right-hand side 
vector for the problem of finding the (k — l)-th derivative. 

So, equation (|40| contains all there is to the sum 5, and the proof is complete. □ 



C Proof of Lemma 3 

We will need the following estimate on S as well 
Repeat of Lemma 3 The sum S — X)j=i satisfies 

I I - Z^l 3 3 I - e n-lf n -k-l)\ 

Proof. Using the formula (JT3J) for Cj and the fact that 



v{5j 



v(S) 



< 



we see that 



(49) 



(50) 



\cj5f\ = k\ 



v{8) 



\o- n -i,n-k-i(S 3 )\\5J\ < k\-—^ 



n — k — 1 



k\A 2n ~ k - 1 fn - l\ A 2n - fe - 1 



e n—l 



(n-l)! 



e™- x (n - fc - 1)! 



and 



11 ™ A2n-t-lf„ i\| A2n-fe-l„| 

V \c-s n \ < V - (n ~ lj! = _ 

I 3 3 I - F n-l( n _ jfe _ I)] e «-l(„ - fc - 1)! 



3=1 



(51) 

(52) 
□ 



D When the Algorithm will Under/Over Estimate 

Sometimes it is useful to know if the algorithm's formula will always underestimate or always overestimate 
the value of f^ k '(t). Here, we show that Theorem [5] gives us the tools to prove the following additional result. 



Repeat of Theorem 3 Assume that all Sj < 

1. If f (n) (t) > 0, then the error term R is positive for all sufficiently small values of h, and the estimate 
of f^ k \t) provided by ([7]) is an underestimate. 



2. If f [n) {t) < 0, then the error term R is negative for all sufficiently small values of h and the estimate 
of f {k) (*) provided by © is an overestimate. 

Proof. Note that if n — k is even, then (— l) n ~ fe = l and <J n ,n-fc(^) > as a sum of products of even 
number of negative S^s. Also, if n — k is odd, then (—l) n ~ h = —1 and o- n , n -k(S) < as a sum of products 
of odd number of negative Si's. In either case, 

(-l) n - k a n ^ k {5) > (53) 

So, if > 0, then 

(-l) n - fc fe"- fc fc!/W(t)g n , n _ fc ffl 

and according to ([20]) . the error term R is positive for all small enough values of h. Similarly, if f (n) (t) < 0, 
then 

{-l) n - k h n - k k\f^{t)a n ^ k (5) ^ Q 
nl 

and according to ([^0)1 , the error term R is negative for all small enough values of h. □ 

E Explicit Form of Coefficients in Special Cases 

We derive some explicit formulas for coefficients produced by algorithm (|13[) in a couple of special cases that 
are common in applications. 

Assume <5j = —i for i = 1, . . . , n, the formula for coefficients Cj simplifies in the case k = 1 to 



and in the case k = 2 to 



c j = (-iy 2 (»)( < ,„, 2 ( 1 ,i,...,I)-I(_i + 1 + i + ... + I] l 



Proof. We start with the general formula 
and observe that for Si = —i the product in the denominator is 



(56) 



= (-j + l)(-i + 2)..-(-l).(l).-.(-j + n) 

= (-ir 1 (j-l)!(n- J )! (57) 



For fe = 1, the symmetric polynomial cr„_x )r j_fc_i becomes 

CT„_l jn _fc_l(5j) = <T n -X tn -2{5j) = —Si ...S n [ - — + — ] (58) 
which with <5; = — z transforms into 



5* ■••*•(-£+£*) 



C r„_ 1 , n _ fc _ 1 (5 i ) = -j(-i) n ^ ( " - I =(-1)"t" ! ( — + >.-) !- r -" = 



Putting the results from ([57)1 , and (jST))) into for fc = 1, we arrive at 
_ (_l)«-i(-l)«i„!(-i+^ =1 i 



(-iy-i(j-l)\(n-j)\ 



(-iy 



j\(n-j)\ y j 



= (-1) J 



which is exactly formula (|54|). 

For fc = 2, the symmetric polynomial cr ra _i n_jt_i becomes 



1 1 



(5 , -i 5 



focusing on the a n -%,2 (ji 



Si ' • ' ' ' <5,_i ' ' • • • ' St. 



i- ) term, we see that 



Cn-1,2 V") 



1 1 



5i £7-1 <5 



CTn,2 











"'0 





which for <5j = —i becomes 
.1 11 

Ofi-1,2 



<5i ' ' <5 7 -i ' S 



= c n _i o 1, • 



1 1 



3 ~ 1 3 + 1 



So putting the results of (|57|) . (jcITj) and ([55]) into (|56p for fc = 2 we obtain 
(-1)"- 2 • 2!i(-l)»n! [<r n , 2 (1, |, . . . , 1) - i (l + • • ■ + 1 - i) 



= (-l)'.2. 7 



(_l)i-i(,--l)!(„-j)! 



3-( n - J)- 

= (-iy-2 



I i 1 1 

2 nj j 



1 1 



II 1 1 
<7 - 2ll '2'---'ny / 



n j 
1 



1 



which is exactly formula ([55 



